We bouwen een populatie met een set aan variabelen waarvan het gemiddelde en spreiding bekend is. Hieruit gaan we dadelijk steekproeven (samples) trekken en onderzoeken welk effect de omvang van de steekproef heeft op de nauwkeurigheid waarmee we een uitspraak kunnen doen over variabele.

Als casus nemen we een ficitieve populatie van 50.000 bankklanten, waarvan we de volgende variabelen weten:

We gaan er vanuit dat dit de hele populatie is (alle klanten van twee banken binnen de onderzochte leeftijdscategorie). In werkelijkheid zal deze populatie onbekend zijn. Organisaties delen dit soort gegevens niet graag.

library(ggplot2)
N <- 50000
sex <- rep(c("M","F"), each=N/2)
age <- rep(18:22, each=N/5)
bank_m <- rep(c("ABN", "RABO"), c(N*0.3,N*0.2))
bank_f <- rep(c("ABN", "RABO"), c(N*0.1,N*0.4))
bank <- c(bank_m, bank_f)
myDF <- data.frame(sex, age, bank)
myDF[sample(1:N,10),]
summary(myDF)
 sex            age       bank      
 F:25000   Min.   :18   ABN :20000  
 M:25000   1st Qu.:19   RABO:30000  
           Median :20               
           Mean   :20               
           3rd Qu.:21               
           Max.   :22               
tbl1 <- table(bank)
prop.tbl1 <- prop.table(tbl1)
p <- prop.tbl1[1]
tbl2 <- table(sex,bank)
tbl2
   bank
sex   ABN  RABO
  F  5000 20000
  M 15000 10000
prop.table(tbl2, 1)
   bank
sex ABN RABO
  F 0.2  0.8
  M 0.6  0.4
prop.table(tbl2, 2)
   bank
sex       ABN      RABO
  F 0.2500000 0.6666667
  M 0.7500000 0.3333333
plot1 <- ggplot(myDF) +
  geom_bar(aes(x=sex, fill=bank))
ggplotly(plot1, width = 800)

plot2 <- ggplot(myDF) +
  geom_bar(aes(x=bank, fill=sex))
ggplotly(plot2, width = 800)

We bouwen nu 3 sets van steekproeven. Iedere set bevat 20 steekproeven. Per set verschilt de steekproefomvang: 10, 50 en 100. Voor iedere steekproef bepalen we de proportie ABN in het totaal. Voor iedere set steekproeven bepalen we de gemiddelde proportie.

generate_samples <- function(N = 1000, n = 10, m = 10) {
  X = matrix(sample(1:N, n*m), n, m)
  return(X)
}
n = 10
m = 20
aselect <- generate_samples(N = N, n = n, m = m)
samples <- sapply(1:m, function(x) myDF[aselect[,x],]$bank)
testABN <- sapply(1:m, function(x) (samples[,x] == "ABN")*1)
propABN <- sapply(1:m, function(x) sum(testABN[,x])/n)
h <- hist(propABN, xlim=range(0,1), col = "mediumaquamarine", main = paste("Histogram of m = " , m, " samples of n = ", n, "items"))
xfit<-seq(min(0),max(1),length=100) 
yfit<-dnorm(xfit,mean=mean(propABN),sd=sd(propABN)) 
yfit <- yfit*diff(h$mids[1:2])*length(propABN) 
lines(xfit, yfit, col="blue", lwd=2)

mean(propABN)
[1] 0.375
sd(propABN)
[1] 0.1585294
(p*(1-p)/n)^0.5
      ABN 
0.1549193 
n = 50
m = 20
aselect <- generate_samples(N = N, n = n, m = m)
samples <- sapply(1:m, function(x) myDF[aselect[,x],]$bank)
testABN <- sapply(1:m, function(x) (samples[,x] == "ABN")*1)
propABN <- sapply(1:m, function(x) sum(testABN[,x])/n)
h <- hist(propABN, xlim=range(0,1), col = "tomato", main = paste("Histogram of m = " , m, " samples of n = ", n, "items"))
xfit<-seq(min(0),max(1),length=100) 
yfit<-dnorm(xfit,mean=mean(propABN),sd=sd(propABN)) 
yfit <- yfit*diff(h$mids[1:2])*length(propABN) 
lines(xfit, yfit, col="blue", lwd=2)

mean(propABN)
[1] 0.415
sd(propABN)
[1] 0.07646809
(p*(1-p)/n)^0.5
       ABN 
0.06928203 
n = 100
m = 20
aselect <- generate_samples(N = N, n = n, m = m)
samples <- sapply(1:m, function(x) myDF[aselect[,x],]$bank)
testABN <- sapply(1:m, function(x) (samples[,x] == "ABN")*1)
propABN <- sapply(1:m, function(x) sum(testABN[,x])/n)
h <- hist(propABN, xlim=range(0,1), col = "tomato", main = paste("Histogram of m = " , m, " samples of n = ", n, "items"))
xfit<-seq(min(0),max(1),length=100) 
yfit<-dnorm(xfit,mean=mean(propABN),sd=sd(propABN)) 
yfit <- yfit*diff(h$mids[1:2])*length(propABN) 
lines(xfit, yfit, col="blue", lwd=2)

mean(propABN)
[1] 0.3905
sd(propABN)
[1] 0.04570788
(p*(1-p)/n)^0.5
       ABN 
0.04898979 

Uit dit experiment blijkt volgende:

Stel we kennen de werkelijke proportie ABN niet. De ABN zelf schat het marktaandeel op p=0.6 We kunnen op basis van de laatste set steekproeven (m=20 en n=100) onderzoeken in hoeverre de schatting realistisch is gegeven de data.

pDist <- 0.6
sdDist <- (pDist*(1-pDist)/n)^0.5
cat("De verdeling van de steekproef gemiddeldes heeft een gemiddelde van", pDist, "en een standaard deviatie van", sdDist, "\n\n")
De verdeling van de steekproef gemiddeldes heeft een gemiddelde van 0.6 en een standaard deviatie van 0.04898979 
pSample <- mean(propABN)
pValue <- pnorm(pSample, mean = pDist, sd = sdDist)

Als de nulhypothese is dat het marktaandeel van de ABN 0.6 bedraagt en we nemen deze voor waar aan dan is de kans dat we een steekproef trekken met een proportie 0.3905 gelijk aan 0.0009%.

Hoe overtuigd ben je na je onderzoek dat de nulhypothese (en de schatting van de ABN) juist is?

LS0tCnRpdGxlOiAiUG9wdWxhdGllcyBlbiBzdGVla3Byb2V2ZW4iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGluY2x1ZGU9RkFMU0UsIHBhZ2VkLnByaW50PUZBTFNFfQpvcHRpb25zKHNjaXBlbj05OTkpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShwbG90bHkpCmBgYAoKV2UgYm91d2VuIGVlbiBwb3B1bGF0aWUgbWV0IGVlbiBzZXQgYWFuIHZhcmlhYmVsZW4gd2FhcnZhbiBoZXQgZ2VtaWRkZWxkZSBlbiBzcHJlaWRpbmcgYmVrZW5kIGlzLiBIaWVydWl0IGdhYW4gd2UgZGFkZWxpamsgc3RlZWtwcm9ldmVuIChzYW1wbGVzKSB0cmVra2VuIGVuIG9uZGVyem9la2VuIHdlbGsgZWZmZWN0IGRlIG9tdmFuZyB2YW4gZGUgc3RlZWtwcm9lZiBoZWVmdCBvcCBkZSBuYXV3a2V1cmlnaGVpZCB3YWFybWVlIHdlIGVlbiB1aXRzcHJhYWsga3VubmVuIGRvZW4gb3ZlciB2YXJpYWJlbGUuCgpBbHMgY2FzdXMgbmVtZW4gd2UgZWVuIGZpY2l0aWV2ZSBwb3B1bGF0aWUgdmFuIDUwLjAwMCBiYW5ra2xhbnRlbiwgd2FhcnZhbiB3ZSBkZSB2b2xnZW5kZSB2YXJpYWJlbGVuIHdldGVuOgoKLSBgc2V4YCAtIGhldCBnZXNsYWNodCB2YW4gZGUga2xhbnQKLSBgYWdlYCAtIGRlIGxlZWZ0aWpkCi0gYGJhbmtgIC0gZGUgbmFhbSB2YW4gZGUgYmFuawoKV2UgZ2FhbiBlciB2YW51aXQgZGF0IGRpdCBkZSBoZWxlIHBvcHVsYXRpZSBpcyAoYWxsZSBrbGFudGVuIHZhbiB0d2VlIGJhbmtlbiBiaW5uZW4gZGUgb25kZXJ6b2NodGUgbGVlZnRpamRzY2F0ZWdvcmllKS4gSW4gd2Vya2VsaWpraGVpZCB6YWwgZGV6ZSBwb3B1bGF0aWUgb25iZWtlbmQgemlqbi4gT3JnYW5pc2F0aWVzIGRlbGVuIGRpdCBzb29ydCBnZWdldmVucyBuaWV0IGdyYWFnLgoKYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQpOIDwtIDUwMDAwCnNleCA8LSByZXAoYygiTSIsIkYiKSwgZWFjaD1OLzIpCmFnZSA8LSByZXAoMTg6MjIsIGVhY2g9Ti81KQpiYW5rX20gPC0gcmVwKGMoIkFCTiIsICJSQUJPIiksIGMoTiowLjMsTiowLjIpKQpiYW5rX2YgPC0gcmVwKGMoIkFCTiIsICJSQUJPIiksIGMoTiowLjEsTiowLjQpKQpiYW5rIDwtIGMoYmFua19tLCBiYW5rX2YpCm15REYgPC0gZGF0YS5mcmFtZShzZXgsIGFnZSwgYmFuaykKbXlERltzYW1wbGUoMTpOLDEwKSxdCnN1bW1hcnkobXlERikKCnRibDEgPC0gdGFibGUoYmFuaykKcHJvcC50YmwxIDwtIHByb3AudGFibGUodGJsMSkKcCA8LSBwcm9wLnRibDFbMV0KCnRibDIgPC0gdGFibGUoc2V4LGJhbmspCnRibDIKcHJvcC50YWJsZSh0YmwyLCAxKQpwcm9wLnRhYmxlKHRibDIsIDIpCgpwbG90MSA8LSBnZ3Bsb3QobXlERikgKwogIGdlb21fYmFyKGFlcyh4PXNleCwgZmlsbD1iYW5rKSkKZ2dwbG90bHkocGxvdDEsIHdpZHRoID0gODAwKQoKcGxvdDIgPC0gZ2dwbG90KG15REYpICsKICBnZW9tX2JhcihhZXMoeD1iYW5rLCBmaWxsPXNleCkpCmdncGxvdGx5KHBsb3QyLCB3aWR0aCA9IDgwMCkKYGBgCgpXZSBib3V3ZW4gbnUgMyBzZXRzIHZhbiBzdGVla3Byb2V2ZW4uIEllZGVyZSBzZXQgYmV2YXQgMjAgc3RlZWtwcm9ldmVuLiBQZXIgc2V0IHZlcnNjaGlsdCBkZSBzdGVla3Byb2Vmb212YW5nOiAxMCwgNTAgZW4gMTAwLiBWb29yIGllZGVyZSBzdGVla3Byb2VmIGJlcGFsZW4gd2UgZGUgcHJvcG9ydGllIEFCTiBpbiBoZXQgdG90YWFsLiBWb29yIGllZGVyZSBzZXQgc3RlZWtwcm9ldmVuIGJlcGFsZW4gd2UgZGUgZ2VtaWRkZWxkZSBwcm9wb3J0aWUuCgpgYGB7cn0KZ2VuZXJhdGVfc2FtcGxlcyA8LSBmdW5jdGlvbihOID0gMTAwMCwgbiA9IDEwLCBtID0gMTApIHsKICBYID0gbWF0cml4KHNhbXBsZSgxOk4sIG4qbSksIG4sIG0pCiAgcmV0dXJuKFgpCn0KCm4gPSAxMAptID0gMjAKYXNlbGVjdCA8LSBnZW5lcmF0ZV9zYW1wbGVzKE4gPSBOLCBuID0gbiwgbSA9IG0pCgpzYW1wbGVzIDwtIHNhcHBseSgxOm0sIGZ1bmN0aW9uKHgpIG15REZbYXNlbGVjdFsseF0sXSRiYW5rKQp0ZXN0QUJOIDwtIHNhcHBseSgxOm0sIGZ1bmN0aW9uKHgpIChzYW1wbGVzWyx4XSA9PSAiQUJOIikqMSkKcHJvcEFCTiA8LSBzYXBwbHkoMTptLCBmdW5jdGlvbih4KSBzdW0odGVzdEFCTlsseF0pL24pCmggPC0gaGlzdChwcm9wQUJOLCB4bGltPXJhbmdlKDAsMSksIGNvbCA9ICJtZWRpdW1hcXVhbWFyaW5lIiwgbWFpbiA9IHBhc3RlKCJIaXN0b2dyYW0gb2YgbSA9ICIgLCBtLCAiIHNhbXBsZXMgb2YgbiA9ICIsIG4sICJpdGVtcyIpKQoKeGZpdDwtc2VxKG1pbigwKSxtYXgoMSksbGVuZ3RoPTEwMCkgCnlmaXQ8LWRub3JtKHhmaXQsbWVhbj1tZWFuKHByb3BBQk4pLHNkPXNkKHByb3BBQk4pKSAKeWZpdCA8LSB5Zml0KmRpZmYoaCRtaWRzWzE6Ml0pKmxlbmd0aChwcm9wQUJOKSAKbGluZXMoeGZpdCwgeWZpdCwgY29sPSJibHVlIiwgbHdkPTIpCgptZWFuKHByb3BBQk4pCnNkKHByb3BBQk4pCihwKigxLXApL24pXjAuNQpgYGAKCmBgYHtyfQpuID0gNTAKbSA9IDIwCmFzZWxlY3QgPC0gZ2VuZXJhdGVfc2FtcGxlcyhOID0gTiwgbiA9IG4sIG0gPSBtKQoKc2FtcGxlcyA8LSBzYXBwbHkoMTptLCBmdW5jdGlvbih4KSBteURGW2FzZWxlY3RbLHhdLF0kYmFuaykKdGVzdEFCTiA8LSBzYXBwbHkoMTptLCBmdW5jdGlvbih4KSAoc2FtcGxlc1sseF0gPT0gIkFCTiIpKjEpCnByb3BBQk4gPC0gc2FwcGx5KDE6bSwgZnVuY3Rpb24oeCkgc3VtKHRlc3RBQk5bLHhdKS9uKQpoIDwtIGhpc3QocHJvcEFCTiwgeGxpbT1yYW5nZSgwLDEpLCBjb2wgPSAidG9tYXRvIiwgbWFpbiA9IHBhc3RlKCJIaXN0b2dyYW0gb2YgbSA9ICIgLCBtLCAiIHNhbXBsZXMgb2YgbiA9ICIsIG4sICJpdGVtcyIpKQoKeGZpdDwtc2VxKG1pbigwKSxtYXgoMSksbGVuZ3RoPTEwMCkgCnlmaXQ8LWRub3JtKHhmaXQsbWVhbj1tZWFuKHByb3BBQk4pLHNkPXNkKHByb3BBQk4pKSAKeWZpdCA8LSB5Zml0KmRpZmYoaCRtaWRzWzE6Ml0pKmxlbmd0aChwcm9wQUJOKSAKbGluZXMoeGZpdCwgeWZpdCwgY29sPSJibHVlIiwgbHdkPTIpCgptZWFuKHByb3BBQk4pCnNkKHByb3BBQk4pCihwKigxLXApL24pXjAuNQpgYGAKCmBgYHtyfQpuID0gMTAwCm0gPSAyMAphc2VsZWN0IDwtIGdlbmVyYXRlX3NhbXBsZXMoTiA9IE4sIG4gPSBuLCBtID0gbSkKCnNhbXBsZXMgPC0gc2FwcGx5KDE6bSwgZnVuY3Rpb24oeCkgbXlERlthc2VsZWN0Wyx4XSxdJGJhbmspCnRlc3RBQk4gPC0gc2FwcGx5KDE6bSwgZnVuY3Rpb24oeCkgKHNhbXBsZXNbLHhdID09ICJBQk4iKSoxKQpwcm9wQUJOIDwtIHNhcHBseSgxOm0sIGZ1bmN0aW9uKHgpIHN1bSh0ZXN0QUJOWyx4XSkvbikKaCA8LSBoaXN0KHByb3BBQk4sIHhsaW09cmFuZ2UoMCwxKSwgY29sID0gInRvbWF0byIsIG1haW4gPSBwYXN0ZSgiSGlzdG9ncmFtIG9mIG0gPSAiICwgbSwgIiBzYW1wbGVzIG9mIG4gPSAiLCBuLCAiaXRlbXMiKSkKCnhmaXQ8LXNlcShtaW4oMCksbWF4KDEpLGxlbmd0aD0xMDApIAp5Zml0PC1kbm9ybSh4Zml0LG1lYW49bWVhbihwcm9wQUJOKSxzZD1zZChwcm9wQUJOKSkgCnlmaXQgPC0geWZpdCpkaWZmKGgkbWlkc1sxOjJdKSpsZW5ndGgocHJvcEFCTikgCmxpbmVzKHhmaXQsIHlmaXQsIGNvbD0iYmx1ZSIsIGx3ZD0yKQoKbWVhbihwcm9wQUJOKQpzZChwcm9wQUJOKQoocCooMS1wKS9uKV4wLjUKYGBgCgpVaXQgZGl0IGV4cGVyaW1lbnQgYmxpamt0IHZvbGdlbmRlOgoKLSBIZXQgZ2VtaWRkZWxkZSAoYG1lYW5gKSB2YW4gZGUgcHJvcG9ydGllIGJpbm5lbiBkZSBhZnpvbmRlcmxpamtlIHNldHMgc3RlZWtwcm9ldmVuIGxpZ3QgZGljaHQgYmlqIGRlIHdlcmtlbGlqa2UgcHJvcG9ydGllOiBgciBwYAotIE5hYXJtYXRlICRuJCB0b2VuZWVtdCwgbmVlbXQgZGUgc3ByZWlkaW5nIHJvbmRvbSBoZXQgZ2VtaWRkZWxkZSBhZjsgZXIga2FuIG5hdXdrZXVyaWdlciB3b3JkZW4gZ2VzY2hhdCB3YXQgZGUgd2Vya2VsaWprZSB3YWFyZGUgdmFuIGhldCBnZW1pZGRlbGRlIGlzLgoKU3RlbCB3ZSBrZW5uZW4gZGUgd2Vya2VsaWprZSBwcm9wb3J0aWUgYEFCTmAgbmlldC4gRGUgQUJOIHplbGYgc2NoYXQgaGV0IG1hcmt0YWFuZGVlbCBvcCBwPTAuNiBXZSBrdW5uZW4gb3AgYmFzaXMgdmFuIGRlIGxhYXRzdGUgc2V0IHN0ZWVrcHJvZXZlbiAobT0yMCBlbiBuPTEwMCkgb25kZXJ6b2VrZW4gaW4gaG9ldmVycmUgZGUgc2NoYXR0aW5nIHJlYWxpc3Rpc2NoIGlzIGdlZ2V2ZW4gZGUgZGF0YS4KCmBgYHtyfQpwRGlzdCA8LSAwLjYKCnNkRGlzdCA8LSAocERpc3QqKDEtcERpc3QpL24pXjAuNQpjYXQoIkRlIHZlcmRlbGluZyB2YW4gZGUgc3RlZWtwcm9lZiBnZW1pZGRlbGRlcyBoZWVmdCBlZW4gZ2VtaWRkZWxkZSB2YW4iLCBwRGlzdCwgImVuIGVlbiBzdGFuZGFhcmQgZGV2aWF0aWUgdmFuIiwgc2REaXN0LCAiXG5cbiIpCgpwU2FtcGxlIDwtIG1lYW4ocHJvcEFCTikKCnBWYWx1ZSA8LSBwbm9ybShwU2FtcGxlLCBtZWFuID0gcERpc3QsIHNkID0gc2REaXN0KQoKYGBgCgpBbHMgZGUgbnVsaHlwb3RoZXNlIGlzIGRhdCBoZXQgbWFya3RhYW5kZWVsIHZhbiBkZSBBQk4gYHIgcERpc3RgIGJlZHJhYWd0IGVuIHdlIG5lbWVuIGRlemUgdm9vciB3YWFyIGFhbiBkYW4gaXMgZGUga2FucyBkYXQgd2UgZWVuIHN0ZWVrcHJvZWYgdHJla2tlbiBtZXQgZWVuIHByb3BvcnRpZSBgciBwU2FtcGxlYCBnZWxpamsgYWFuIGByIHJvdW5kKHBWYWx1ZSwgNikqMTAwYCUuCgoqKkhvZSBvdmVydHVpZ2QgYmVuIGplIG5hIGplIG9uZGVyem9layBkYXQgZGUgbnVsaHlwb3RoZXNlIChlbiBkZSBzY2hhdHRpbmcgdmFuIGRlIEFCTikganVpc3QgaXM/KioK